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We search for di-muon decays of a low-mass Higgs boson (^4°) produced in radiative T(1S) decays. 
The T(1S) sample is selected by tagging the pion pair in the Y(2S, 3S) — >• ir + ir~Y(lS) transitions, 
using a data sample of 92.8 x 10 6 T(2S) and 116.8 x 10 6 T(3S) events collected by the BABAR 
detector. We find no evidence for A production and set 90% confidence level upper limits on the 
product branching fraction B(T(1S) -> 7 A ) x B{A° ->• in the range of (0.28 - 9.7) x 10 -6 

for 0.212 < m A o < 9.20 GeV/c 2 . The results are combined with our previous measurements of 
T(2S, 3S) —¥ jA°, A —¥ fi + \i~ to set limits on the effective coupling of the 6-quark to the A . 
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Many extensions of the Standard Model (SM), such 
as the Next-to-Minimal Supersymmetric Standard Model 
(NMSSM), include a light Higgs boson [IS]. The Mini- 
mal Supersymmetric Standard Model (MSSM) 3] solves 
the hierarchy problem of the SM, but fails to explain why 
the value of the /i-parameter is of the order of the elec- 
troweak scale, which is many orders of magnitude below 
the next natural scale, the Planck scale. The NMSSM 
solves this so-called "^-problem" [3] by adding a singlet 
chiral superfield to the MSSM, generating an effective jU- 
tcrm. As a result, the NMSSM Higgs sector contains a 
total of three neutral CP-even, two neutral CP-odd, and 
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two charged Higgs bosons. The lightest CP-odd Higgs 
boson (A ) could have a mass smaller than twice the 
mass of the 6-quark [l] , making it detectable via radia- 
tive r(nS) -> jA° (n = 1,2,3) decays @. 

The coupling of the ^4° field to up-type (down-type) 
fermion pairs is proportional to cos^a cot/3 (cos#Atan/3), 
where #a is the mixing angle between the singlet com- 
ponent and the MSSM component of the A , and tan/3 
is the ratio of the vacuum expectation values of the up 
and down type Higgs doublets. The branching fraction of 
Y(IS) —> jA° could be as large as 10~ 4 depending on the 
values of the ^4° mass, tan/? and cos6*a [2j. Constraints 
on the low-mass NMSSM Higgs sector are also important 
for interpreting the SM Higgs sector 

BaBAR has previously searched for ^4° production 
in several final states |7Hig|. including T(2S, 35) — > 
■yA°,A Q —> fi + n~ Q. Similar searches have been per- 
formed by CLEO in the di-muon and di-tau final states 
in radiative T(1S) decays [ll|, and more recently by BE- 
SIII in J/V> j A , A -> fi+fj,- [13, and by the CMS 
experiment in pp — > A , A — > /i + /i~ 13]. These results 
have ruled out a substantial fraction of the NMSSM pa- 
rameter space 
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We report herein a search for a di-muon resonance 
in the fully reconstructed decay chain of T(2S, 35) — > 
7r+7r-r(15), T(15) -> 7 A°, A ->■ /i+pt-. This search 
is based on a sample of (92.8 ± 0.8) x 10 6 T(25) and 
(116.8 ± 1.0) x 10 6 T(35) mesons collected by the BABAR 
detector at the PEP-II asymmetric-energy e + e _ collider 
located at the SLAC National Accelerator Laboratory. 
A sample of T(15) mesons is selected by tagging the 
di-pion transition, which results in a substantial back- 
ground reduction compared to direct searches of A in 
T(25,35) -> 7A decays. We assume that the light 
Higgs boson that we search for is a scalar or pseudoscalar 
particle with a negligible decay width compared to the 
experimental resolution [Tij . 

The BABAR detector is described in detail elsewhere 
[H, Q3] • Charged particle momenta are measured in a 
five-layer double-sided silicon vertex tracker and a 40- 
layer drift chamber, both operating in a 1.5 T solenoidal 
magnetic field. Charged particle identification (PID) is 
performed using a ring-imaging Cherenkov detector and 
the energy loss (dE/dx) in the tracking system. Photon 
and electron energies are measured in a CsI(Tl) electro- 
magnetic calorimeter, while muons arc identified in the 
instrumented magnetic flux return of the magnet. 

Monte Carlo (MC) simulated events are used to study 
the detector acceptance and to optimize the event selec- 
tion procedure. The EvtGen package [l8| is used to simu- 
late the e + e~ — > qq (q = u, d, s, c) and generic T(2S, 35) 
production, BHWIDE to simulate the Bhabha scat- 
tering, and KK2F [2(| to simulate the processes e + e~ — > 
and e + e~ — > (-f)r + T~. Dedicated MC sam- 
ples of T(2S, 35) generic decays to 7r + 7r~T(15) with 
T(15) — s- 7/i + ^i~ decays, hereafter refered to as the 
"non-resonant di-muon decays" are also generated. Sig- 
nal events are generated using a phase-space (P-wave) 
model for the A -> /Lt+^t" (T(15) -> jA°) decay, and the 
hadronic matrix elements measured by the CLEO experi- 
ment H3 are used for the T(25, 35) ->■ 7r+7r-r(15) mod- 
eling. The detector response is simulated by GEANT4 
[22! ] and time-dependent detector effects are included in 
the simulation. A sample corresponding to about 5% of 
the dataset is used to validate the selection and fitting- 
procedure. To avoid bias, this sample is discarded from 
the final dataset. We perform a blind analysis, where 
the rest of the T(25, 35) datasets are blinded until the 
analysis procedure is frozen. 

We select events containing exactly four charged tracks 
and a single energetic photon with a center-of-mass (CM) 
energy larger than 200 MeV. The tracks are required to 
have a distance of closest approach to the interaction 
point of less than 1.5 cm in the plane transverse to the 
beam axis and less than 10 cm along the beam axis. At 
least one of the tracks must be identified as a muon by 
particle ID algorithms; the probability for misidentify- 
ing a charged pion as a muon is 3%. Additional photons 
with CM energies below the threshold of 200 MeV are 
also allowed to be present in the events. The two highest 
momentum tracks in the CM frame with opposite charge 



are assumed to be muon candidates and are required to 
originate from a common vertex to form the ^4° candi- 
dates. 

The T(15) candidate is reconstructed by combining 
the A candidate with the energetic photon candidate 
and by requiring the invariant mass to be between 9.0 
and 9.8GeV/c 2 . The T{2S, 35) candidates are formed 
by combining the T(15) candidate with the two re- 
maining tracks, assumed to be pions. The di-pion in- 
variant mass must be in the range [2m„ (jnT(2S,3S) ~ 
mrns))], compatible with the kinematic boundaries of 
the T(25,35) -> tt + tt-T{1S) decay. The entire decay 
chain is then fit by imposing the decay vertex of the 
T(25, 35) candidate to be constrained to the beam inter- 
action region, and a mass constraint on the T(1S) and 
T(25, 35) candidates, as well as requiring the energy of 
the T(25, 35) candidate to be consistent with the e + e~ 
CM energy. These constraints improve the resolution of 
the di-muon invariant mass to less than 10 MeV/c 2 . 

To improve the purity of the T(15) sample, we train 
a Random Forest (RF) classifier [23| on simulated signal 
and background events, using variables that distinguish 
signal from background in the T(2S, 35) — > ir + ir~T(lS) 
transitions. The following quantities are used as inputs 
to the classifier: the cosine of the angle between the two 
pion candidates in the laboratory frame, the transverse 
momentum of the di-pion system in the laboratory frame, 
the azimuthal angle and transverse momentum of each 
pion, the di-pion invariant mass (m n7V ), the pion helicity 
angle, the transverse position of the di-pion vertex and 
the mass recoiling against the di-pion system, defined as 

nirccoii = V s + m l-n - 2v^ES^> where y/s is the e+e" 
CM energy and E„„ is the CM energy of the di-pion 
system. For signal-like events, the m reC oii distribution 
peaks at the mass of the T(15), with a mass resolution of 
about 3 MeV/c 2 . The RF output peaks at 1 for signal- like 
candidates and peaks at for the background-like can- 
didates. The optimum value of the RF selection is cho- 
sen to maximize Punzi's figure-of-merit, e/(0.5N a + yB) 
[2H , where N a = 3 is the number of standard deviations 
desired from the result, and e and B are the average ef- 
ficiency and background yield over a broad m^o range 
(0.212 - 9.20 GeV/c 2 ), respectively. 

A total of 11,136 T(25) and 3,857 T(35) candidates 
are selected by these criteria. Figures [T] and [2] show the 
distributions of the m re coil and di-muon reduced mass, 
Tired = \J ra 'l l + t _ l - ~ 4m 2 , together with the background 

prediction estimated from the MC samples, which are 
dominated by the non-resonant di-muon decays. The 
reduced mass is equal to twice the momentum of the 
di-muon system in the rest frame of the ^4°, and has a 
smooth distribution in the region of the kinematic thresh- 
old TUu+a- ~ 2m^ (r7i rc d w 0). After unblinding the 
data, two peaking components corresponding to p° and 
J/tp mesons are observed in the T(3S) dataset. The 
/9°-mesons are mainly produced in initial state radiation 
(ISR) events, along with two or more pions. This peak 
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FIG. 1: (color online) The distribution of m re coii for (a) 
the T(2S) and (b) the Y(3S) datasets, together with the 
background predictions from the Monte Carlo (MC) samples, 
which are dominated by the non-resonant di-muon decays. 
The MC samples are normalized to the data luminosity. 



disappears if we require both candidates to be identi- 
fied as muons in the A reconstruction. The J/ip mesons 
arise from e + e _ — > jjsriP{2S), tp(2S) — s- tt + tt~ J/ip , 
J/i/j — > p + p~ decays. The J/ip and p° events in the 
T(2S) dataset are suppressed since the di-pion mass dis- 
tribution in these events is above the kinematic edge of 
the di-pion mass distribution of T(2S) decays, but well 
within the range of values allowed for the Y(3S) decays. 

We extract the signal yield as a function of m A o in 
the region 0.212 < m A o < 9.20 GeV/c 2 by perform- 
ing a series of one-dimensional unbinncd extended max- 
imum likelihood (ML) fits to the m red distribution. We 
fit over fixed intervals in the low mass region: 0.002 < 
TO rcd < 1-85 GeV/c 2 for 0.212 < m A o < 1.50 GeV/c 2 , 
1.4 < m rod < 5.6 GeV/c 2 for 1.50 < m A o < 5.36 GeV/c 2 
and 5.25 < m rcd < 7.3 GeV/c 2 for 5.36 < m A a < 7.10 
GeV/c 2 . Above this range, we use sliding intervals 
pL - 0.2 GeV/c 2 < m rod < Li + 0.15 GeV/c 2 , where pL 
is the mean of the reduced mass distribution. 

The probability density function (PDF) of the signal 
is described by a sum of two Crystal Ball (CB) func- 
tions [25j. The signal PDF is determined as a function 
of m A o using signal MC samples generated at 26 differ- 
ent masses, and by interpolating the PDF parameters be- 
tween each mass point. The resolution of the m rGd distri- 
bution for signal MC increases monotonically with m A o 
from 2 to 9 MeV/c 2 , while the signal efficiency decreases 
from 38.3% (40.4%) to 31.7% (31.6%) for T(2S) (T(3S)) 



FIG. 2: (color online) The distribution of m rc d for (a) the 
T(2S) and (b) the T(3S) datasets, together with the back- 
ground predictions from the various Monte Carlo (MC) sam- 
ples. The MC samples are normalized to the data luminosity. 
Two peaking components corresponding to the p° and J/tp 
mesons are observed in the Y(3S) dataset. The contribution 
of J/tjj mesons is not included in the MC predictions, whereas 
the p° meson is poorly modeled in the MC. 



transitions. The background for m^o < 1.5 GeV/c 2 is 
described by a threshold function 

/(m rcd ) cx [Erf(s(m rod -m )) + l]+cxp(c + cim rod ) (1) 

where s is a threshold parameter and mo is determined 
by the kinematic end point of the m rod distribution, and 
Co and ci are the coefficients of the polynomial function. 
These parameters are determined from the MC sample of 
the non-resonant di-muon decays. In the range of 1.5 < 
771^0 < 7.1 GeV/c 2 , the background is modeled with a 
second order Chebyshev polynomial function and with a 
first order Chebyshev polynomial function for m^o > 7.1 
GeV/c 2 . We model the p° background with a Gaussian 
function, using the sideband data of the di-pion recoil 
mass distribution from the T(3S) sample to determine 
its mean and width. The background in the J/tp mass 
region is modeled by a CB function using a MC sample of 
e+e- j IS r^(2S), ip(2S) ir+n-J/ip, J/tp p+ p~ 
decays. 

We search for the A signal in steps of half the m rod res- 
olution, resulting in a total of 4,585 points. The shape of 
the signal PDF is fixed while the continuum background 
PDF shape, the signal and background yields are allowed 
to float. In the fits to the T(3S) dataset, we include a 
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J/tp and p° background components whose shapes are 
fixed, but we allow their yields to float. To address the 
problem associated with the fit involving low statistics, 
we impose a lower bound on the signal yield by requir- 
ing that the total signal plus background PDF remains 
non- negative [Hj]. 

The results of the fits for some selected m A a points 
are shown in Fig. [3j Figure [4] shows the number of 
signal events (N s i g ) and the signal significance (S) as 
a function of m A o. The signal significance is defined 
as <S = sign(N sig )^~ 2 ln(£ /£ max) i where £>rnax is the 
maximum likelihood value for a fit with a floating signal 
yield centered at m A o , and Cq is the likelihood value for 
N S ig — 0. The largest values of significance are found to 
be 3.62ct at m A o = 7.85 GeV/c 2 for T(2S) dataset, 2.97cr 
at m A o = 3.78 GeV/c 2 for Y(3S) dataset and 3.24cr at 
m A a = 3.88 GeV/c 2 for the combined Y(2S, 3S) dataset. 
We estimate the probability of observing a fluctuation of 
S > 3.62cr (S > 2.97ct) in the Y{2S) (T(35)) dataset to 
be 18.1% (66.2%), and for S > 3.24ct in the combined 
T(2S, 3S) dataset to be 46.5% based on a large ensemble 
of pseudo-experiments. Therefore, the distribution of the 
signal significance is compatible with the null hypothesis. 

Tables|I]and [TI]summarize the additive and multiplica- 
tive systematic uncertainties, respectively considered in 
this analysis. Additive uncertainties arise from the choice 
of fixed PDF shapes and from a possible bias on the fit- 
ted signal yield. The systematic uncertainty associated 
with the fixed parameters of the PDFs is determined by 
varying each parameter within its statistical uncertainties 
while taking correlations between the parameters into ac- 
count. This uncertainty is found to be small for each 
mass point and does not scale with the signal yields. We 
perform a study of a possible fit bias on the signal yield 
with a large number of pseudo-experiments. The biases 
are consistent with zero and their average uncertainty is 
taken as a systematic uncertainty. 

The multiplicative systematic uncertainties arise from 
the signal selection and T(nS) counting. They in- 
clude contributions from the RF classifier selection, par- 
ticle identification, photon selection, tracking and the 
Y(2S, 3S) constrained fit. The uncertainty associated 
with the RF classifier is studied using the non-resonant 
di-muon decays in both data and MC. We apply the RF 
selection to these control samples to calculate the rela- 
tive difference in efficiency between data and MC. The 
systematic uncertainty related to the photon selection is 
measured using an e + e~ — > 77 sample in which one of 
the photon converts into an e + e~ pair in the detector 
material [9j. The uncertainties on the branching frac- 
tions Y{2S, 3S) -> 7r+7r-T(lS')_are 2.2% and 2.3% for 
the T(2S) and Y(3S) datasets [27|, respectively. 

We find no significant signal and set 90% confidence 
level (C.L.) Bayesian upper limits on the product of 
branching fractions of B(Y(1S) -> -fA°)xB(A° -> ^ + M~) 
in the range of 0.212 < m A a < 9.20 GeV/c 2 . Figure [5] 
shows the branching fraction upper limits at 90% C.L. 
which are determined with flat priors. The systematic 




m red (GeV/c0 



FIG. 3: (color online) Projection plots from the Id unbinned 
ML fit to the m re d distribution in the Y(2S) dataset for 
m A o = 5.93 GeV/c 2 (top) and m A o = 7.85 GeV/c 2 (bot- 
tom) that return the largest upward fluctuations. The green 
dotted line shows the contribution of the signal PDF, the ma- 
genta dashed line shows the contribution of the continuum 
background PDF and solid blue line shows the total PDF. 
The signal peak corresponds to a statistical significance of 
3.38<r for the top plot and 3.62<r for the bottom plot. Based 
on the trial factor study, we interpret these results as mere 
background fluctuations. 



uncertainty is included by convolving the likelihood with 
a Gaussian distribution having a width equal to the sys- 
tematic uncertainties described above. The combined re- 
sult is obtained by simply adding the logarithms of the 
Y(2S, 35) likelihoods. The limits range between (0.37 — 
8.97) x 10~ 6 for the Y(2S) dataset, (1.13 - 24.2) x 10~ 6 
for the Y(3S) dataset, and (0.28 - 9.7) x 10~ 6 for the 
combined Y(2S, 3S) dataset. 

The branching fractions of B(Y(nS) —> jA°) (n = 
1,2,3) are related to the effective Yukawa coupling (fx) 
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FIG. 4: (color online) The number of signal events {N s i g ) 
and significance obtained from the fit as a function of m^o 
for (a,b) the T(2S) and (c,d) the T(3S) datasets. The impact 
of the requirement that the total PDF remains non-negative 
is clearly visible in the lower m A o region. 



of the &-quark to the A" through d [H 

B{T(nS) -> 7 A°) _ fj f m\ \ 
B(r(nS) -> 1+1-) 2na\ m 2 ,^/' [) 

where I = e or fj, and a is the fine structure constant. The 
value of fr incorporates the QCD and relativistic correc- 
tions to B{T(nS) — > yA ) [29l|. as well as the leptonic 
width of T(nS) — > l + l~ [30] - These corrections are as 
large as 30% to first order in the strong coupling constant 
(as), but have comparable uncertainties |31j . The 90% 
C.L. upper limits on /£ x B(A° —> fj, + n~) for combined 
Y(2S, 3S) datasets range from 0.54 x 10~ 6 to 3.0 x 10" 4 
depending upon the mass of the A , which is shown in 



TABLE I: Additive systematic uncertainties and their 
sources. 

Source Y(2S) (events) T(3S) (events) 

Nsi g PDF (0.00 - 0.62) (0.04 - 0.58) 

Fit Bias 0.22 0.17 

Total (0.22 - 0.66) (0.17 - 0.60) 



TABLE II: Multiplicative systematic uncertainties and their 
sources. 



Source 


r(2S) (%) 


r(35) (%) 


Muon-ID 


4.30 


4.25 


Charged tracks 


3.73 


3.50 


B{T{nS) -> 7r+7r-T(lS)) 


2.20 


2.30 


RF classifier 


2.21 


2.16 


Photon efficiency 


1.96 


1.96 


T(nS) kinematic fit \ 2 


1.52 


2.96 




0.86 


0.86 


Total 


7.00 


7.32 



Fig. Eljd). We combine these results with our previous 
measurements of T(2S,3S) -> jA°, A -> 0, to 

obtain 90% C.L. upper limits on j\ x B(A° ->■ in 
the range of (0.29 - 40) x 10~ 6 for m A o < 9.2 GeV/c 2 
(Fig. He)). 

In summary, we find no evidence for a light scalar 
Higgs boson in the radiative decays of T(1S) and set 
90% C.L. upper limits on the product branching frac- 
tion of B{T{IS) -> jA°) x B(A° -> n+n~) in the 
range of (0.28 - 9.7) x 10" 6 for 0.212 < m A a < 9.20 
GeV/c 2 . These results improve the current best limits 
by a factor of 2-3 for m^o < 1.2 GeV/c 2 and are com- 
parable to the previous BaBAR results [2fl in the mass 
range of 1.20 < m A a < 3.6 GeV/c 2 . We also set lim- 
its on the product x B(A° — > /i + /i _ ) at the level of 
(0.29 - 40) x 10~ 6 for m A o < 9.2 GeV/c 2 . 
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FIG. 5: (color online) The 90% C.L. upper limits (UL) on the 
product of branching fractions B(T(1S) -> -yA°) x B{A° -> 
fi + H~) for (a) the T(2S) dataset, (b) the T(3S) dataset and 
(c) the combined Y(2S, 3S) dataset; (d) the 90% C.L. UL on 
the effective Yukawa coupling f r X B(A° fi + together 
with our previous BABAR measurement of T(2S, 3S) — > 7A , 
A" —¥ )i + 0], and (e) the combined limit. Details of the 
UL and Yukawa coupling as a function of m A o are provided 
in II. 
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